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Abstract 

We present a new isogeometric method for the discretization of the Reissner- 
Mindlin plate bending problem. The proposed scheme follows a recent theo- 
retical framework that makes possible to construct a space of smooth discrete 
deflections Wh and a space of smooth discrete rotations 0^ such that the 
Kirchhoff contstraint is exactly satisfied at the limit. Therefore we obtain a 
formulation which is natural from the theoretical/mechanical viewpoint and 
locking free by construction. We prove that the method is uniformly stable 
and satisfies optimal convergence estimates. Finally, the theoretical results 
are fully supported by numerical tests. 

Keywords: Isogeometric analysis, Reissner Mindlin plates, De Rham 
diagram 



1. Introduction 

The Reissner-Mindlin theory is widely used to describe the bending be- 
havior of an elastic plate loaded by a transverse force. Despite its simple 
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formulation, the discretization by means of finite elements is not straightfor- 
ward, since standard low-order schemes exhibit a severe lack of convergence 
whenever the thickness is too small with respect to the other characteris- 
tic dimensions of the plate. This undesirable phenomenon, known as shear 
locking, is nowadays well understood: as the plate thickness tends to zero, 
the Reissner-Mindlin model enforces the Kirchhoff constraint, which is typ- 
ically too severe for Finite Element Methods (FEM), especially if low-order 
polynomials are employed (see, for instance, the monograph by Brezzi and 
Fortin Roughly speaking, the root of the shear locking phenomenon is 
that the space of discrete functions which satisfy the Kirchhoff constraint is 
very small, and does not properly approximate a generic plate solution. The 
most popular way to overcome the shear locking phenomenon in FEM is to 
reduce the influence of the shear energy by considering a mixed formulation 
and/or suitable shear reduction operator. As a consequence, the choice of 
the discrete spaces requires particular care, also because of the possible oc- 
currence of spurious modes. A vast engineering and mathematical literature 
is devoted to the design and analysis of plate elements for FEM; we mention 
here, in a totally non-exhaustive way, the works 0, IE 0, 0, @, Q, 0, 0, 0, II 
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In 2005, IsoGeometric Analysis (IGA) has been introduced by T.J.R. 
Hughes and co-authors in [l3| as a novel technique for the discretization of 
partial differential equations. IGA is having a growing impact on several 



fields, from fluid dynamics 14|, [15[, to structural mechanics [16|, [17|, llSl, [19 
and electromagnetics 0, 21] . A comprehensive reference for IGA is the book 

IGA methodologies are designed with the aim of improving the inter- 
operability between numerical simulation of physical phenomena and the 
Computer Aided Design (CAD) systems. Indeed, the ultimate goal is to 
drastically reduce the error in the representation of the computational do- 
main and the re-meshing by the use of the "exact" CAD geometry directly 
at the coarsest level of discretization. This is achieved by using B-Splines 
or Non Uniform Rational B-Splines (NURBS) for the geometry description 
as well as for the representation of the unknown fields. The use of Spline 
or NURBS functions, together with isoparametric concepts, results in an ex- 
tremely successful idea and paves the way to many new numerical schemes 
enjoying features that would be extremely hard to achieve within a standard 
FEM. Splines and NURBS offer a flexible set of basis functions for which 
refinement, de- refinement, degree elevation and mesh deformation are very 
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efficient (e.g., [HI). Beside the fact that one can directly treat geometries 
described by Splines and NURBS parametrizations, these functions are in- 
teresting in themselves since they easily allow global smoothness beyond the 
classical C°-continuity of FEM. This feature has been advantageously ex- 
ploited in recent works, for example 17 1, and [23[, and studied in 24 1, 25 



Furthermore, the intrinsic regularity of the Spline basis functions opened 
the way to completely new discretization schemes for Maxwell equations 



20l.l21|. as well as for other problems, such as, the Stokes problem [15(. These 
schemes are based on suitable smooth approximation of differential forms, 
verifying a De Rham diagram in the spirit of [26]. The proposed theoretical 
framework can be usefully adopted also for discretizing the Reissner-Mindlin 
plate system which is the object of the present paper. Indeed, it makes 
possible to construct a space of smooth discrete deflections Wh as discrete 
0-forms and a space of smooth discrete rotations 0^ as discrete 1-forms such 
that it holds: 

VW h c & h . 

This allows us to select approximation spaces containing a subspace which 
both exactly satisfies the Kirchhoff constraint and has optimal approximation 
features. Therefore, we obtain a simple formulation which is natural from 
the theoretical/mechanical viewpoint since it is built to be locking free, and 
it is easy to study. 

Finally, the method inherits all the advantages, mentioned above, which 
are typical of IGA: the capability to incorporate exactly CAD geometries, 
and the flexibility in the choice of the polynomial degree and function reg- 
ularity. Regarding the first advantage, the importance of reproducing the 



exact geometry in plate analysis has been underlined very clearly in [27f . 

The outline of the paper is the following. Section [2] is devoted to the de- 
scription of Reissner-Mindlin model and E] provides basics definition of spline 
spaces. Our discretization method is proposed in Section H] and analysed in 
Section [51 Finally, Section [6] is devoted to the numerical validation. 



2. The Reissner-Mindlin plate bending problem 

Let Q be a Lipschitz bounded open set of R 2 representing the midsurface 
of the plate and let T be its boundary. Other assumptions on the domain Q 
will be set at the end of Section 13.21 We assume that the boundary T is the 
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union of three disjoint sets 



r = r / ur s ur c , 

with r/,r s ,r c being a finite union of connected components. The plate 
is clamped in r c , simply supported in T s and free in Tf. We assume for 
simplicity of exposition that all boundary conditions are homogeneous, and 
that the union T s U T c has positive measure, in order to neglect rigid body 
motions of the plate. Finally, let T s be divided into a soft and a hard part 
(both being a finite union of connected components), i.e. T s = T ss UT^. We 
then set the spaces for deflections and rotations 

W = {v G H\n) : v = 0on T S UT C } 

& = {rie[H 1 (n)} 2 : r/ = on r c , T] ■ t = on T sh }, 

with t the unit tangent to T obtained by an anti-clockwise rotation of the 
outward normal n. 

Following the Reissner-Mindlin model, see for instance jlj], the plate bend- 
ing problem requires to solve 



Find e O , w eW such that 

a(0, rj) + fikt~ 2 (0 - Vw, rj-Vv) = (f, v) V77 G ©o, v G W 



where /1 is the shear modulus and k is the so-called shear correction factor. 
Above, t represents the plate thickness, w the deflection, 6 the rotation of 
the normal fibers and / the applied scaled transversal load. Moreover, (-, ■) 
stands for the standard scalar product in L 2 (Q) and the bilinear form a(-, ■) 
is defined by 

a(0,T7) = (C £ (0), £ (77)), 

with C the positive definite tensor of bending moduli and e(-) the symmetric 
gradient operator. Introducing the scaled shear stresses 7 = fikt~ 2 (6 — 
Vw), Problem ([T]) can be written in terms of the following mixed variational 
formulation: 



' Find 6 G ©0, w G W , 7 G S such that 
a(e,ri) + (rr,ri-Vv) = (f,v) \/rie@ ,veW 



t 2 



(2) 



;e-Vro,s)--(7,s) = VsGS, 
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where £ = [L 2 ({1)} 2 . 

The above problem is well posed and t-uniformly stable when using the 
H 1 norm for the spaces ©o, Wo, and the norm: 

, , || | (s,77 — Vv) 

s s := ^||s||l 2 + sup 



(r,,v)e® xw (WvWm + IMI^i) 1/2 ' 

for the space £ (see [lj], for instance). To simplify notation, and without any 
loss of generality, we will assume //& = 1 in the analysis that follows. 

3. B-splines and piece-wise smooth functions 

3.1. B-spline spaces and piece-wise smooth functions in one dimension 

Given positive integers p and n, such that n > p + 1, we introduce the 
ordered knot vector 

£ := {0 = £l> £2, • • • , £n+p+l = 1} ) 

where we allow repetitions of knots, that is, we assume £1 < £2 < ■ ■ ■ < 
£ n+p+ i. In the following we will only work with open knot vectors, which 
means that the first p + 1 knots in 5 are equal to 0, and the last p+1 
are equal to 1, and we assume that all internal knots have multiplicity r, 
1 < r < p + 1, so that 




; Cl > C2, • - ; C2 , • • • , Cm, • ■ ■ , Cm) - 



r times p+1 times 

The vector 

-2 = {0 = Ci,Ca,...,Cm = l} 

represents the (ordered) vector of knots without repetitions, and the relation 
m = + 2 holds. 



Through the iterative procedure detailed in [13] we construct p-degree 
(that is, (p + l)-order) B-spline basis functions, denoted by B i: for i = 
1, . . . , n. These basis functions are piecewise polynomials of degree p on the 
subdivision {£1, . . . , ( m }. At Q they have a := p — r continuous derivatives. 
Therefore, — 1 < a < p — 1: the maximum multiplicity allowed, r —p + 1, 
gives a = — 1, which stands for a discontinuity at each Q. Each basis function 
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Bi is non- negative and supported in the interval [£i,£i+ p +i]- Moreover, these 
B-spline functions constitute a partition of unity, that is 

n 

^Bi(x) = l VxG(0, 1). 

i=l 

The space of B-splines spanned by the basis functions will be denoted by 

SI := spa n {^}f =1 . 

Derivatives of splines are splines as well. Let S^Xi an d be spline spaces 
constructed, according with the notation above, on the same subdivision 
{Ci> • • • ^ Cm}- Then, it is easy to see that 

Notice moreover that 

#^=p+l + (m-2)(p-Qt), (4) 

and 

#S p a + + \=p + 2 + (m-2)(p-a), (5) 

where # is used to denote the dimension of the linear space. Then, from 
O^©, ^S'a+i = ft^a + 1> m agreement with the fact that the derivative is 
a surjective operator from to and has a one-dimensional kernel, the 
constants. 

We denote by the space of piecewise smooth functions on (0, 1), whose 
restriction to each subinterval (Q, admits a C°° extension to the closed 
interval [d, Ci+i] with a continuous derivatives at Q, for all i = 2, . . . , m — 1. 

3.2. B-spline spaces and piece-wise smooth functions in two dimensions 

The definition of B-splines spaces given above can be extended to two di- 
mensions as follows. Let us consider the square Q = (0, l) 2 C M 2 , which will 
be referred to as parametric domain. Given integers p d , r d , n d and a d = p d — 
r d , with d = 1, 2, we introduce the knot vectors E d = {£ M , £ 2 ,d, £n d + Pd +i,d} 
and the associated vectors Z d = {Ci,<z> ■ ■ • , Cm d ,<i} as i n the one- dimensional 
case. Associated with these knot vectors there is a mesh of the parametric 
domain, that is, a partition of (0, l) 2 into rectangles: 

&h = {Q = ®d=i,2(Ci d ,d, d d +i,d), ^<id<m d - 1}. (6) 
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Given an element Q G Qh, we set Hq = diam(Q), and define the global mesh 
size h = max{/iQ, Q G f^}. 

In this paper we make the following assumption: 

Assumption 3.1. The parametric mesh Qh is shape regular. 

It is important to remark that, due to the tensor product structure, 
shape regularity implies quasi-uniformity, i.e., there exists a positive 
constant k, fixed once and for all, such that 

kh<h Q <h , VQ G n h . (7) 

We associate to the two given knot vectors E d , d = 1,2 the p^-degree 
univariate B-splines basis functions B^, with i = 1, . . . , n^. Then, on the 
associated mesh Qh, we define the tensor-product B-spline basis functions as 

Bij := B iA <g> B jt2 , i = 1, ■ ■ ■ , n h j = 1, . . . , n 2 . 

Then, the tensor product B-spline space is defined as the space spanned by 
these basis functions, namely 

S p a \% = S£S(n fc ) := S£ ® S "l = Bpan^}^ . 

Notice that the space <S**'^ i^h) is fully characterized by the mesh f^> by 
Pi, P2, oci and a 2 , as our notation reflects. The minimum regularity of the 
space is a := min{ad : d — 1, 2}. 

In a similar way, we define on the space of piecewise smooth functions 
with interelement regularity on the vertical and horizontal mesh edges given 
by «! and a 2 respectively. This is denoted by 

r°° - r°° (O, ) — C°° a C°° 

Precisely, a function in a2 admits a C°° extension in the closure of each 
element Q G i\, has oti continuous derivatives on the edges {(xi,x 2 ) : a?i = 
Ci,i, 0,2 < £2 < (7+1,2}, for j = 1, . . . , m 2 - 1, i = 2, . . . , m 1 - 1 and a 2 
continuous derivatives on the edges {(xi,x 2 ) : (7,1 < x i < Cj+i,i' x 2 — 0,2, }, 
for j = 1, . . . , mi — 1, i = 2, . . . , m 2 — 1. From the definitions, C a2 . 

From an initial coarse mesh f2^ , refinements are constructed by knot in- 
sertion (with possible repetition, see j2f|). Therefore, we end up considering 
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a family of meshes {Qh}h<h an d associated spaces, with the global mesh size 
h playing the role of family index, as usual in finite element literature. 

We assume that our computational domain Q C M 2 can be exactly 
parametrized by a geometrical mapping F : Q — > Q which belongs to 
(C™ a2 (flh)) 2 , with piecewise smooth inverse, and is independent of the mesh 
family index h. The geometrical map F naturally induces a mesh Qh on 
Q, which is the image of ([6]). Notice that, as a conseguence of ([7]) and the 
boundness of F and its inverse, there exist two positive constants k', k", 
fixed once and for all, such that 

k'h < h K < k"h VK G Qh ■ (8) 

Finally , Tf, T s T c C T = <9f2 denote the preimage of Tf T s and T c , 
respectively. 

4. Discretization of the scalar and vector fields 

4-1- Spline spaces on the parametric domain 

Given the two (horizontal and vertical) knot vectors, Z\ = {Ci,i, • • • , Cmi,i} 
and Z 2 = {Ci,2, • • • , (m 2 ,2}, and the associated mesh Q h on the parametric do- 
main Q, with 1 < a < p + 1, we introduce the following spaces: 

W h = S?%(h h ) (9a) 

® h = s p a --Z& h ) x s%i{n h ) (9b) 

S h = VW h + @ h (9c) 

Since by construction VWh Q &h, it clearly holds that = 0^. 

Essential boundary conditions has to be set directly in the function spaces 
and this brings us to the definition: 

W h>0 = W h , (p, a) = {ve S p a %(Q h ) : v = on f s U f c } (10a) 

@h,o = 0h,o(P 5 ") = {^7 e S p a Z\%(n h ) x S p a p a Z\(h h ) : r; = on f c and T] ■ t = on f s/l } 

(10b) 

S^o = S fc>0 (p, a) = VW hfl + 0^,0 (10c) 

where i is a unit tangent vector at T sh . Notice that now VWh,o Q ®h,o 
only if r c = and T ss = (see e.g., [20]). As a conseguence, H ht0 needs a 
characterization which is the object of the next Lemma: 
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Lemma 4.1. The following characterization holds 
se @ h 

s G Sfc o <3> { s-t = OonT sh l>T c 



s(x) = 0, for every corner x G dQ such that x G T ss fl T c . 

_ (") 
Moreover, given s G S^ci; ^ ere are 11 Wa,o> *7 e ®M> s,uc ^ 

77 — Vf = s, 

IMIif^n) + IMIh 1 ^) - ^ H s llL 2 (n)- 

Proof. Denote S^o the space of fields that fulfill the characterization on the 
right hand side of (ITT]) . Let S^q 3 s = 77 s — Vf s . Then, from ([3]), we easily 
see that Vv s belongs to S^Z\' p a (^h) X ^'(O/,), and the boundary condi- 
tions that define fllOaj) — fllOcD implies that r/ s — Vv s fulfills the homogeneous 
boundary conditions stated in (ITTI) . Thus E^o C E^o- 

In order to prove the inclusion E^o ^ ^^,0; notice that the fields in E^o 
which vanish on r c belong to @h,o, therefore we need to show that the B- 
spline basis functions of S^ q whose support intersect T c belong to VW^o- 
This is again easily seen from the structure of the spaces. 

Let again E^o 3 s = rj s — Vv s . Define v G Wh,o such that 

= —s ■ n on T c (13) 

and such that the B-spline coefficients of v not involved in ( I13p are set to 
zero. Because of the locality of the B-spline basis functions, v is supported 
in a strip of element associated with the p + 1 knot spans around T c . By 
construction, the norm of v is bounded by the norm of ( ITS"]) , then, by a scaling 
argument and inverse estimate, 



= Ch 1 ' 2 

<c\ 



dv 



d n L 2 (r c ) 

II ( 14 ) 

S ' U llL2(f c ) 



lL 2 (H) • 

Since s + Vt> vanishes on T c , we can set 

ry = s + Vv G Gk o, 
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and, using again an inverse estimate and (j!4p . the following estimate holds 

II^IIh 1 ^) - ^ IMIl 2 (6) 

= Ch- 1 (||s|| i2( n) + l|Vu|| z2(fi) J (15) 
< Ch 1 \\s\\ L 2(ny 

The property (|T2l) follows from (|T4|) and (|T5l) . □ 

^.2. Spline spaces on the physical domain 

Once the finite dimensional spaces Wh^, @h,o and S^ q on the parametric 
domain Q have been defined, we construct the corresponding spaces Wh t o, 
@h,o and S^o in the physical domain Q. 

The deflections space Wh,o is mapped from the reference domain via the 
geometrical parametrization F : Q — > Q, that is 

W hfi = {voF~ 1 :v e W h}0 (p,a)}. 

For the rotations and shears spaces we use the covariant map: 

M = {DF- T v o F^ 1 : v G M }, 

S M = {DF- T v o F- 1 : v G 

where DF~ T = (DF _1 ) T is the transpose of the gradient of the inverse 
of the geometrical map F . The push-forward by the covariant map has 
two main properties: i) it preserves the nullity of tangential components, 
ii) it maps gradient to gradients, i.e., by the chain rule, we have Vu> = 
D F " T Vw o F _1 when w = w o F -1 . Thus, we have that v ■ t = on 
for all v G ®h,o and, recalling Lemma I4.1[ the following holds: 

Lemma 4.2. We have 

£h,o = W M + M , 
and, for s G T, h , there exist v G Wh 0; rj G ©/, 0; such that 



r\ — Vf = s, 

IvWrnin) + Ikllif^n) < C^ _1 ||s||l 2 (^)- 



(16) 
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We now turn to approximation properties and we make use of the ap- 
proximation results proved in 29] and adapted in jiof for vector fields under 



covariant transformation. 

Lemma 4.3. Let 1 < a < p + I, and let F G {C^ a {Vt h )f . There exist 
projectors U Wh : W ->■ W hfi , Tl &h0 : 6 ->■ fti o and n Sfe 
that, for all 1 < s < p + 1 

- n WM >)||tf 1(K) < c/^M^, wk e n h ,wv ew n h s (k), 

for all 1 < s < p 

h - n eM (»7)||Hi(K) < Ch'- x \-n\ H . {k) , G O h , Vry G O n (# S (K)) 2 

for all 1 < s < p 

\\ S - n SM ( s )|| i2(K) < c^M^, vat g Vs g s s , 

E s = (& + VW )n(H s (K)) 2 . 
Proof. The projector Hw h can be chosen as the quasi-interpolant proposed 



m 



291 ] and homogeneuous boundary condition are imposed by setting to zero 
che coefficient of all the basis functions which are non-zero on T c and T s . The 
projector N-& h0 is the one constructed in [20], with boundary condition set 



only on T c and T sh ; and IT Sh is the same with the set of boundary condition 
corresponding to ffTTj) . □ 

Remark 4.1. // T c = and T ss = 0, then VW hfi C @ hfi = T, hfi . In this 
case, we can choose as interpolation operators the commuting projectors that 
are constructed in Jg^J: i.e., there are projectors P\v h0 and P@ h0 such that 



We refer to ]M, \2QJ for the details 



4-3. The discrete problem 

We can now introduce our proposed method. The Galerkin formulation 
for the approximation of Problem (CQ), reads 

{Find (9 h ,Wh) G ®h,o x W fti0 such that 
a(0 h , r} h ) + r 2 (0 h - Vw h , rj h - Vv h ) = (f, v h ) \/(r} h , v h ) G Q h<0 x W hfi ; 

(17) 
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this is the discrete problem which is solved in our isogeometric code. For the 
purpose of its stability and convergence analysis, we introduce the discrete 
shear stress 

lh = r\e h - Ww h ) e s M , 

and reformulate Problem ([IT]) in the following mixed form (analogous to 
Problem (ED): 



Find d h w h ,~f h ) G M X W h)0 x £ M such that 
a{d h ,r) h ) + ( lh1 ri h -Vv h ) = (f,v h ) Vri h G @ h ,o,v h G W hfi (18) 
(0 h - \7w h , s h ) - t 2 (7 h , s h ) = Vs ft G T, hfi . 

Introducing, for all (/3, u, r) and (77, v, s) in © x Wo x S, the symmetric 
bilinear form 

B((3, u, r; 17, u, a) = a(/3, 77) + (r, 77 - Vv) + (/3 - Vw, s) - t 2 (r, a), (19) 
the continuous and discrete problems can be also written as 

B(d J w J r,r),v,s) = (f,v) V(77,^,a)G0 o x^ o xS, (20) 

and 

B(d h , w h , -y h ; rj h , v h , s h ) = (/, u h ) V^, u/j, a fc ) G M x W hfi x £ M . (21) 

Note that, differently from what happens with most finite element meth- 
)ds (e.g. the well-known MITC elements, cf. 30, 5[), if one considers the limit 



oc 

case t — in formulation (fT8"j) . the second equation gives exactly 6 h = V%. 
Therefore, substituting such identity in the first equation of (fIBl) . one gets 
the problem 

Find Wh G Wh,o such that 
a(Vw h , Vv h ) = (/, G W^ )0 , 

i.e. the Kirchhoff plate bending problem. 

5. Convergence analysis 

In this section we investigate the convergence properties of the method 
proposed in Section HI 
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Under the assumption T c = T ss = 0, we present in Section loTTl an analysis 
based only upon coercivity. A similar approach is proposed by Duran and 
Libermann in Q, but things are straightforard in our case, thanks to the 
commuting projectors that are at our disposal (see Remark 14. ip . 

On the other hand, for general boundary conditions, we need a more 
delicate analysis which is the object of Section 15.21 

5.1. Coercivity-based convergence analysis 
We assume in this section: 

Assumption 5.1. It holds T c = and T ss = 0. 

We recall that, under the above assumption, as noted in Remark 14.11 the 
following properties hold 

©ft,0 = £fc,o , -P® h , V = VP Wh0 . 

The following lemma holds. 

Lemma 5.1. Let (0,w) be the solution of Problem ([T]) and (0h,Wh) be the 
solution of Problem ffTTp . Let the shears 7 = t~ 2 (0 — Vw) and *y h = t~ 2 (0h — 
Vwh) ■ Then for 1 < s < p it holds 

\\0-0h\\m(n)+\\w-w h \\ H i in) +t\\j-j h \\ L 2 {n) < C/i^^H^H^^+tllTll^- 1 ^)) , 

(22) 

where the constant C is independent of h,t. 
Proof. From ([1]) and (jTTl) . we obtain 

a(0 - 6 h , r) h ) + (7 - 7 h , r) h - Vv h ) = \/r) h e & hl v h G W h . 

Therefore, for every (0i,wi) G @o,h X Wo,h, and jj := t~ 2 (0i — Vwi) G 
So,h, we infer 



a(0j - h , -q h ) + (7 7 - 7 h , r) h - Vv h ) = a(0! - 0, rj h ) 

+ hi - 7, Vh ~ Vv h ) Vrj h e@ h ,v h eW h . 

Choosing r\ h = 0j — 0h and Vh = wi — Wh, we get 
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t 2 {li-lh) =ri h -Vv h . 
Therefore, from (|23|) we get 

a {e t - e h , 0/ - e h ) + t 2 ( 7/ - 7fc , 7j - 7fc ) = o (0j - 0, 0, - e h ) 

+ * 2 (7/-7,7/-7/ 1 )- 

(24) 

Using standard arguments, from (|24p we obtain 

\\9i - dh\\m{n) + t\hi - lh\\L 2 (n) < C(\\6 I - 0\\ H i^ + £|| 7j - tI U 2 (n)) • 

We choose now 0/ = P® h0 6 and u>/ = PyK h0 u '> where -P© h0 and Pw h0 w are 
the interpolation operators introduced in Section 14.21 By the choice of 7/ , 
we obtain: 

7/ := t-\0j - Vwj) = r 2 P &hfi (6 - Vw) = P &hfil . 
Therefore, also using Lemma [4. 3[ we have 

*ll7/ - 7lU 2 (n) = *||-Pe kl0 7 ~ 7lU 2 (n) < Cti , ~H\\y\\ H —i(ny (25) 



Furthermore, again due to Lemma [4. 3 [ it holds 

||0/ — ^Wh 1 ^) = \\P® h>0 ~ 0|U 1 (^) < Ch s ||0||H«-i(n). (26) 
Using the triangle inequality and (f2"5l) - (l26l) . from Lemma I5TT1 we get 

). (27) 

To get the error estimate on the deflections, we simply notice that it holds 

V(w - w h ) = (0 - h ) - t 2 ( 7 - lh ). (28) 
Estimate ((22J now easily follows from (|27|)-(|2g|). 
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5.2. Convergence in a discrete norm 

In this section we suppose that general boundary conditions are set and 
thus Assumption 15.11 does not hold. On the other hand, we remind that (jSJ) 
holds true and this is crucial in the subsequent analysis. 

We will make use of the following discrete norm j^]: 

Willie = Nlfri(n)+ Yl t 2 + fe 2 W^v - r)\\ 2 L2{K) 



Kcn h ~ ' K 



t'\\s\\' r ^n\-\- ^2 ^k\\ S \\l 2 (K) (}J ' 



|2 _ +2|| c ||2 
Is — 1 \\ S \\L 2 (Q.) 



Ken 



h 



iii^^^iii 2 = mollis + ni s iiis > 

for all r] G ©o, v G Wq and sGS. It is easy to check that it holds 
1 1 1 7 ?) ^Ille > C\\v\\ 2 Hl{n) Vry G @o, v G Wq , 

with C independent of h and t. Therefore, the norm ||| • \ also includes 
the H 1 norm of both rotations and deflections. 

The following stability result can be shown adopting well-known tech- 
niques, cf. jsj, combined with the results Section HI Therefore, the proof will 
be only sketched and several details will be omitted. 

Proposition 5.1. For all ((3 h , u h , r h ) G Q ,h, W ,h, ^o,h, there exists (rj h , v h , s h )) G 
®o,h, W ,h, S ,h such that 

B((3 h ,u h ,T h ;ri h ,v h ,s h ) > d\ \ \(3 h , u h , r h \ | | 2 , (30) 
\\\r] h ,v h ,s h \\ \ < C 2 \\\f3 h ,u h ,T h \\\ , (31) 

with C\ and positive constants independent of h and t. 

Proof. Taking v\, s\) = ((3 h , — r^) in ffl9|) . one immediately gets 

B({3 h ,u h ,T h ;r)l,vlsl)>C\\p h \\ 2 Hl(n) +t 2 \\T h \\ 2 L 2 (n) , (32) 
111*71^4111 < \\\Ph,Uh,T h \\\ . (33) 
Due to ( |T6l) . we can choose (rj^v 2 , s\) G ©o,/i x Wq^ X Hh,o such that 

4 = o 

rf h - Vv 2 h = h 2 r h (34) 
H^hlli/ 1 ^) + Ikhllff 1 ^) — Ch\\r h\\L 2 (n) ■ 
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The Cauchy-Schwarz inequality yields 

B{f3 h ,u h ,T h ;ril,v 2 h ,s 2 h ) = a{f3 h ,rf h ) + h 2 \ \r h \ || 2(n) 

> h2 \\ T h\\ 2 L 2(n) - C\\^ h \\ H i{a)\\n\\\m{^) (35) 

> h 2 \\r h\\ 2 L 2 (n) ~ Ch\\(3 h \\ H i^)\\T 2 h \\ H i {n ) . 

Applying the arithmetic-geometric mean inequality, from (I34p - (l35p one easily 

gets 

B((3 h ,u h} T h -r 1 2 h y h ,s 2 h ) > -h 2 \\T h \\ 2 L2{n) - C'\\f3 h \\ 2 Hl(n) . (36) 

Furthermore, by definition (1291) . recalling ([7]) and one has 

( h 2 \ 

1 1 1^,^X1 1 1 < C[h\\T\\ L 2 (n) + — _||T||£8 ( n)J < C\\\(3 h ,u h ,T h \\\ . (37) 

We finally select {rf h ,vl, s\) = (0, 0, (t + h)- 2 (/3 h - Vu h )). Using again the 
arithmetic-geometric mean inequality and some simple algebra, we obtain 

1 t 2 
B{P h ,u h ,T h -rf h ,vl,sl) = +h y \\Ph- Vu h \\ L 2 m - y_ + h y (T,P h -Vu h ) 

1 t 4 



1 

2(t 2 + h 2 ) 



> nU9. , u9J Ph-Vy>h\\ 2 L H n)-Ct 2 \\T h \\ 2 Lm . 



(38) 

Moreover, recalling (129]) and ((7|), it follows 

\\\vlvls 3 h \\\ ^Cit + h^WPh-VuhWvw <C\\\(3 h ,u h ,T h \\\ . (39) 
We now consider the linear combination 

3 



{r} h ,v h ,s h ) = J^c^X^ 



with positive q G M. Estimates f l3"Uj) and f l3"Tj) follow from a suitable choice 
of the coefficients c h by using (152"])-(133Th (|MD (EZD and (J3ED-(EHD- 

n 



From Proposition I5.1[ the following result is easily deduced. 
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Proposition 5.2. Let (0,w) be the solution of Problem (pQ) and (0h,Wh) be 
the solution of Problem (1171) . Let the shears 7 = t~ 2 (0 — Vw) and 'y h = 
t~ 2 (dh — Vwt). Then, for any (0j,wi,jj) in ©o,/i x Wo,h x So,h it holds 

\\\0 - h ,w - w h ,~/ - 7JH < C\\\6 - 0i, w - u;/,7 - 7j||| , 

where the constant C is independent of h,t. 

Proof. The proof is a consequence of standard consistency-stability argu- 
ments. We apply Proposition [5J] to the difference (0h — 0i,Wh — wi,j h — ^j), 
then apply the error equation given by (1201) and (T2TT) . to derive 

1 1 10^-0/, Wh-wi, 7/1-7/1 1 1 2 < CB(0-O I ,w-w I ,j-j I ;r} h ,v h ,s h ) , (40) 

where 

|||»7h,Vft, s/illl < Hl^fc - - w/>7h-7illl • ( 41 ) 

Observing that the bilinear form B is bounded with respect to the norm 1 1 1 • 1 1 
and applying (f4"T|) . bound (HOT) gives 

1110^ - 0j,w/i - io/, 7 ft - 7/1 1 1 < C|||0 - 0/,w - if/, 7- 7/| 1 1 . 
The result follows using the triangle inequality. 

□ 

Combining the above proposition with the interpolation results of Lemma 14731 
the next Corollary easily follows. 

Corollary 5.1. Let (0,w) be the solution of Problem (|T]) and (0h,Wh) be 
the solution of Problem ( 1171) . Let the shears 7 = t~ 2 (0 — Vw) and *y h = 
t~ 2 (0h — Vwh)- Then, provided the solution of (pQ) is sufficiently regular for 
the right-hand-side to make sense, for all 1 < s < p it holds 

1110-0^,^-^,7-7^111 < C/l^ 1 ^||0|| s + ||w|| s+ l + ||7l|max(O, S -2)+t||7l|s-l) , 

where the constant C is independent of h and t. 

We note that the regularity required in Corollary 15.11 may be unrealistic 
for high values s, due to the presence of boundary layers (see for instance 

ft- 
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5.3. Improved convergence rates 

We start this section with the following lemma, which shows an improved 
convergence property for the L? norm of the rotations. We omit the proof, 
since it merely adopts a classical Aubin-Nietsche argument. 

Lemma 5.2. Let the Problem (CQ) be regular. Then, under the notation of 
Proposition ^. 2l it holds 

\\0 - 6h\\L 2 (n) < Ch{\\0 - 0fc||ffi(n) + i||7 — 7/JU 2 (fi)J 
with C independent of h,t. 

We remark that the regularity of Problem (CQ) holds, for instance, when 
the domain Q is convex and the plate is clamped on the whole boundary (32| . 
For more general cases we refer the reader to [ill, 32] where the (smooth) 
boundary layers and the effect of corners are respectively investigated. 

The following result gives an improved convergence rate of convergence 
for the deflection variable. 

Lemma 5.3. Under the assumptions and notation of Lemma \5.°A for any 
interpolant wi G Wh,o, it holds 

\\w-w h \\ H i(n) < Ch(\\0 - 0fc||jri(n) + *H7-7fclU»(n)J + ||V(u> - w 7 )|| L 2 (n) 
with C independent of h and t. 

Proof. Testing Problems ([TJ and (fT7|) with the choice (0, Vh) G Q^xW^q, 
and taking the difference, one obtains 

(V(tw - w h ), Vv h ) = {6- 9 h , Vv h ) . 

Thus, setting Vh = uih — Wi one gets 

\\V(w h - w/)||ia ( n) = (V(wfc-ioi),Vv/,) 

= (V(w h - w), Vv h ) + (V(w - iuj), Vv h ) (42) 
= (0 - h , Vv h ) + (V(tw - wj), Vv A ). 



Applying the Poincare inequality, then using (|42p and the Cauchy-Schwarz 
inequality, we obtain 

\\wh-wi\\w(n) < C\\V(w h -wj)\\ L 2 {Q) < C\\0-0 h \\ L 2 {Q) + \\V(w-wj)\\ L 2 {n) . 
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The result trivially follows from the above bound, combined with the triangle 
inequality and Lemma 15.21 

n 

Note that in the proposed method the approximation order of the space 
Wh is one point higher with respect to the rotation space 0^, as shown in 
Lemma 14.31 Therefore, the above result indeed implies an improved rate of 
convergence with respect to Proposition 15.21 



6. Numerical tests 

In this section we present some numerical experiments to show the actual 
performance of the numerical methods introduced in Section |H For all tests, 
we select a material with Poisson's ratio v = 0.3 and Young's modulus E = 
1.092 • 10 7 N/m 2 . Accordingly, we will always express forces and lengths in 
N and m, respectively. 

In all the tests we start from the mesh used for the geometry represen- 
tation and we build the different components of the approximation spaces 
performing three simple steps: 

1. degree elevation of the basis function used for the mapping, to the 
requested degree of the given field component; 

2. knot repetitions in order to have the proper continuity; 

3. the h- refinement is performed using uniform knot insertions in the orig- 
inal knots vectors. 

We remark that, if the geometry is described by NURBS, step 1. is 
perfomed on the B-splines basis obtained by setting the NURBS weight to 
one. 

6.1. Case 1: fully hardly clamped square plate. 



We consider a problem having a known the analytical solution (see 33j). 
This test consists of a unitary square block [0, l] 2 with all the four sides 
clamped, i.e., Y = T c , subject to a body load given by 

f{x, y) = 12{ f_ y2) [I2y(j/ - l)(5x 2 - 5x + l)(2y 2 (y - l) 2 

+ x(x - l)(5y 2 - 5y + 1)) 

+ 12x{x - l)(5y 2 -5y + l)(2x 2 (x - l) 2 

+ y(y - l)(5x 2 -5x + 1))] . 
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The analytical solution is 



0{x,y) 



V(z/-l)V(x-l) 2 (2x-l) N 
x\x-lfy\y-l) 2 {2y-l) j 



w(x,y) = ^s 3 (x-l)V(2/-l) 3 
2t 2 

- — v [y 3 (y - l f< x - W 5 * 2 - 5x + 1) 

5(1 — v) 

+ x\x - l) 3 y(y - l)(5y 2 - 5y + 1)] . 

In Figure[T]we plot the if 1 -norm approximation error for the displacement 
field w(x, y) and for the rotation field 6(x, y). We run tests with p = 3, 4 ad 
a = 2, 3 both for deflections and rotations, respectively, and the thickness 
t = 10~ 3 . Figure [T] clearly displays convergence rates in perfect agreement 
with the theoretical results of Section [51 

In addition, we performed simulations (not reported here) using also 
t G {1CT 1 , 1CT 2 , 1CT 4 }, and we found that the methods are substantially 
insensitive to thickness changes. 

6.2. Case 2: quarter of an annulus with four hardly simply supported sides. 

The second test consists in a fully hardly simply supported quarter of an 
annulus plate., i.e., Y = T sh (see Figure [2]) and loaded by 

f(x,y) = 10 4 sin (2 arctan(|)) . 

In this case the analytical solution is not available. Therefore, we use as 
reference solution the one computed on a fine mesh (h = 1/256) and with 
degree p = 3 and regularity index a = 2 both for deflections and rotations. 
The same degree and regularity index are used for the simulations on the 
coarser meshes. 

In Figure [3] and Figure H] we plot the if 1 -norm approximation errors for 
t = 10 -2 and t = 10~ 3 respectively. Once again, an excellent agreement with 
the theoretical predictions can be noticed. 
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Figure 1: Case 1 (square with four sides hardly clamped), t = 10 




0.5 1.0 1.5 2.0 2.5 

Figure 2: Fully hardly simply supported annular plate. 
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Figure 3: Case 2 (fully hardly simply supported annular plate), p = 3, t = 10 2 . 

6.3. Case 3: boundary layer test problem 

The third test problem shares the same load and geometry used for Case 2 
(Figure [2]), but with different boundary conditions. More precisely, we set 
r c = and 

• T ss = {(x,y) ER 2 : x 2 + y 2 = l}; 

• T sh = {(0, y) : 1< y < 2.5} U {(x, 0) : 1< x < 2.5} 
. r, = {(a;,y)GM a : x 2 + 2/ 2 = f}. 

On the curved parts of the boundary, this problem exhibits a boundary layer 



whose characteristic length is 0(t) (see 3l|). In our computations we set 



t = lCT 2 , and p = 3, a = 2 both for deflections and rotations. For this test 

1 /2 

case we always plot, in log-log scale, the errors versus N^ OF , where N DO f 
denots the total number of degrees of freedom. We remark that for uniform 
meshes, N^q F behaves like hr 1 . 
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As for Case 2, the analytical solution is not available, and a reference 
numerical solution is obtained using a very fine mesh. 

Figure E] displays the error behaviour when using a uniform mesh refine- 
ment. It is worth noticing that a severe suboptimal convergence rate occurs 
as long as the mesh is not sufficiently fine to resolve the boundary layer. Of 
course, this phenomenon implies that uniform meshes requires an excessive 
number of degrees of freedom to significantly reduce the error, especially for 
small plate thicknesses. 

We now employ a sequence of meshes adapted to the boundary layers 
(Figure E]), using the following procedure. An initial mesh consisting of three 
elements is set up as in Figure [5j The width of both the layer elements is 
0.045 (corresponding to 3% of the total annulus width). The finer meshes 
are then obtained by uniform refinement of each of the three initial elements. 
Figure [7] shows the error behaviour using that mesh sequence. We notice 
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0.5 1.0 1.5 2.0 2.5 

Figure 5: Coarsest mesh for the layers-adapted case. 




Figure 6: Case 3 (quarter of ring with boundary layers). Uniform mesh, p = 3, t = 10 
that the optimal convergence rate for the if 1 -norm error is now restored. 
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Figure 7: Case 3 (quarter of ring with boundary layers). Layers-adapted mesh, p = 3, 

t = io- 2 . 
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